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HIGHLIGHTS 


•  A  2-d  transient  mathematical  model  to  simulate  the  electrochemical  lithiation  in  silicon  nano-wire  structure  is  developed. 

•  The  model  geometry  is  an  axisymmetric  cylindrical  silicon  nanowire  electrode  anchored  to  a  copper  current  collector  substrate. 

•  Model  physics  include  diffusive  transport  of  lithium,  volume  expansion,  and  diffusion  induced  stresses. 

•  Larger  radius  of  the  nanowire,  and  faster  lithiation  rates  lead  to  increased  radial  and  tangential  stresses  in  the  nanowire. 
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This  paper  presents  a  2-dimensional  transient  numerical  model  to  simulate  the  electrochemical  lithium 
insertion  in  a  silicon  nanowire  (Si  NW)  electrode.  The  model  geometry  is  a  cylindrical  Si  NW  electrode 
anchored  to  a  copper  current  collector  (Cu  CC)  substrate.  The  model  solves  for  diffusion  of  lithium  in  Si 
NW,  stress  generation  in  the  Si  NW  due  to  chemical  and  elastic  strains,  stress  generation  in  the  Cu  CC  due 
to  elastic  strain,  and  volume  expansion  in  the  Si  NW  and  Cu  CC  geometries.  The  evolution  of  stress 
components,  i.e.,  radial,  axial  and  tangential  stresses  in  different  regions  in  the  Si  NW  are  presented  and 
discussed.  The  effect  of  radius  of  Si  NW  and  lithiation  rate,  on  the  maximum  stresses  developed  in  the  Si 
NW  are  also  discussed. 

©  2014  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Silicon  electrode  is  pursued  as  a  potential  negative  electrode  for 
lithium-ion  batteries  owing  to  its  high  gravimetric  (mAh  g-1)  and 
volumetric  capacity  (mAh  L-1)  compared  to  the  existing  state  of  the 
art  graphite  electrode  [1].  One  of  the  critical  challenges  in  the 
commercialization  of  Si  electrode  is  to  minimize  particle  fracture 
developed  during  lithiation  and  delithiation  of  the  Si  electrode 
[2,3].  Recent  experimental  studies  have  demonstrated  the  use  of 
nano-size  Si  structures  as  electrodes.  These  electrode  structures 
exhibited  minimal  particle  fracture  and  also  enabled  repeated 
cycling  [4-6].  While  different  mechanisms  have  been  proposed  for 
this  behavior,  a  detailed  physics  based  analysis  combining  the 
electrochemical  and  structural  aspects  of  lithium  insertion  in  such 
nanostructures  have  not  been  undertaken.  A  variety  of  detailed 
phenomenological  models  exists  in  the  literature  for  lithium 
intercalation  in  porous  electrodes,  which  treat  the  transport  of 
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electrolyte  due  to  diffusion  and  migration,  reaction  kinetics  at  in¬ 
terfaces,  and  transport  of  Li  and  electrons  in  solid  phase  [7-13].  The 
general  modeling  framework  presented  in  these  papers  cannot  be 
directly  used  to  simulate  advanced  high  capacity  electrodes,  spe¬ 
cifically  the  alloy  type  electrodes  such  as  Si,  Sn  etc.,  because  (a)  the 
stresses  developed  during  lithium  insertion/deinsertion  and  (b) 
volume  change  associated  with  lithium  insertion/deinsertion  are 
not  considered.  So  to  accurately  model  such  type  of  high  capacity 
electrodes  which  undergo  substantial  volume  changes,  particle 
level  expansion/contraction  and  electrode  level  displacement  along 
with  build  up  of  stresses  have  to  be  captured.  Early  research  by 
Prussin  [14]  demonstrated  that  the  diffusion  induced  stresses 
generated  by  concentration  distributions  are  of  similar  nature  to 
the  thermal  stresses  developed  in  an  elastic  medium.  Modeling  of 
diffusion  induced  stresses  was  also  studied  in  detail  by  other  re¬ 
searchers  [15-18]  for  different  geometries  such  as  hollow  cylinders, 
plates  etc.  Similar  approaches  were  extended  to  battery  electrode 
chemistries  on  a  particle  level  to  calculate  intercalation  induced 
stresses  assuming  no  volume  changes.  Zhang  et  al.  [19]  presented  a 
numerical  model  to  calculate  diffusion  induced  stresses  for  spher¬ 
ical  and  ellipsoidal  shaped  LiM^CU  single  particle.  Also,  the  work 
by  Cheng  and  Verbrugge  [20,21]  derived  analytical  expressions 
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(assuming  negligible  pressure  induced  diffusion  and  no  volume 
change)  to  calculate  stresses  that  arise  from  concentration  gradi¬ 
ents  for  a  spherical  particle.  This  modeling  framework  was  also 
incorporated  into  a  porous  electrode  framework  [22].  All  the  above 
referenced  models  in  addition  to  other  published  work  [23-25] 
were  developed  assuming  dilute  solution  theory  for  diffusion 
within  particle,  with  no  moving  boundaries  (negligible  volume 
change)  and  for  low  expansion  materials. 

Christensen  et  al.  [26]  presented  a  more  rigorous  mathematical 
framework  based  on  concentrated  solution  theory,  which  in¬ 
corporates  volume  expansion  and  stress  build  up  in  a  single 
spherical  particle  electrode  and  case  studies  for  lithiation  in  a 
spherical  carbon  particle  (8%  volume  expansion)  were  discussed. 
The  same  framework  was  also  used  to  calculate  the  stresses  in 
LiMn204  single  spherical  particle  electrode  27]  and  was  also  later 
extended  to  porous  electrodes  [28]  containing  graphitic 
mesophase-carbon-microbead  (MCMB)  anode  and  lithium  man¬ 
ganese  oxide  spinel  cathode.  The  author  also  emphasizes  the 
importance  of  thermodynamic  factor,  pressure  driven  diffusion  and 
extent  of  volume  change  in  determining  the  cell  potential  profiles 
and  initiation  of  fracture.  However  most  of  the  above  mentioned 
work  was  based  on  electrodes  which  undergo  volume  change  in  the 
order  of  10%.  To  model  large  volume  expansion  in  electrodes, 
Chandrasekaran  et  al.  [29]  modeled  a  single  particle  Si  electrode 
under  galvanostatic  and  potentiodynamic  control  of  lithiation  of  Si 
to  Li3.75Si  associated  with  a  270%  volume  change.  In  a  later  paper 
[30],  the  same  approach  was  extended  to  a  porous  electrode  to 
describe  how  particle  level  expansion  affects  the  porosity  of  the 
electrode.  The  authors  ignored  stress  calculations  based  on  the 
assumption  that  the  nano  sized  particles  would  not  build  appre¬ 
ciable  concentration  gradients  to  generate  diffusion  induced 
stresses.  Gao  et  al.  [31  ]  modeled  stress  build  up  due  to  concentra¬ 
tion  gradients  for  a  1-d  (radial)  cylindrical  geometry  for  a  nano 
sized  Si  electrode  for  a  dilute  solid  solution  with  constant  density. 
The  authors  also  discuss  the  strong  coupling  between  stress 
enhanced  diffusion  and  diffusion  induced  stresses  for  electrodes 
associated  with  large  volume  expansion. 

In  this  paper  (Part  I),  we  present  a  model  to  describe  diffusion 
and  stress  build  up  in  a  2-d  silicon  nano  wire  (Si  NW)  geometry 


anchored  to  a  Cu  substrate  under  galvanostatic  conditions.  The 
model  in  general  follows  the  framework  described  in  reference  [26] 
but  applied  to  Si  electrode  with  a  maximum  lithiation  to  Li3.7sSi 
associated  with  a  270%  volume  change.  Case  studies  related  to 
lithiation  rate  &  nanowire  radius  are  presented.  In  the  second  part 
of  the  paper  (Part  II),  a  dimensionally  reduced  1-d  model  (in  the 
radial  co-ordinate)  is  presented  and  the  predictions  between  the 
1-d  and  2-d  model  are  compared.  The  concentration  and  stress 
profiles  predicted  by  the  1-d  model  compares  well  with  the  2- 
d  model  at  regions  far  away  from  the  Cu  substrate.  Consequently, 
the  1-d  model  was  used  to  analyze  different  geometries  -  nanotube 
&  core/shell  structures  and  the  evolution  of  stresses  in  these  ge¬ 
ometries  are  compared  to  the  nanowire  geometry. 

2.  Model  development 

2.1.  Assumptions 

The  geometry  of  the  Si  NW  anchored  to  the  Cu  CC  substrate  is 
shown  in  Fig.  1  (left);  the  initial  unexpanded  radius  and  the  length 
of  the  Si  NW  are  Ruw  =  50  nm  and  HNW  =  10  p  respectively.  The 
model  geometry  consists  of  a  2-d  axisymmetric  cut  from  the  overall 
geometry  as  shown  in  Fig.  1  (right),  wherein  the  dependence  of  the 
variables  in  the  6  direction  was  ignored.  Other  key  assumptions  in 
the  model  were 

•  The  charge  storage  mechanism  in  the  Si  electrode  is  modeled  by 
considering  the  electrochemical  reaction  of  Li  at  the  surface  of 
the  Si  NW  followed  by  transport  of  Li  into  the  Si  NW.  The 
lithiated  Si  mixture  is  represented  as  a  solid  solution,  therefore 
phase  transitions  are  ignored. 

•  The  transport  of  the  host  silicon  is  solely  due  to  the  convective  flux, 
where  as  the  transport  of  lithium  is  due  to  the  combined  effect  of 
gradients  in  concentration  and  pressure,  and  convective  flux. 

•  Lithium  diffusion  into  Cu  CC  substrate  is  ignored  and  therefore 
the  Si  NW/Cu  CC  interface  acts  as  a  Li  blocking  interface. 

•  Stress  strain  relationship  was  assumed  to  obey  Hooke’s  law 
(linear)  in  the  entire  lithiation  regime.  For  particles  of  nano¬ 
meter  radii,  the  stress  generated  due  to  insertion  is  typically  less 


Fig.  1.  Schematic  of  the  Si  NW  anchored  to  the  Cu  CC  substrate  (left).  2-d  axisymmetric  slice  of  Si  NW  anchored  to  Cu  CC  base  substrate  used  as  the  geometry  for  the  2-d  model 
(right). 
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than  the  yield  stress  limit  for  onset  of  plasticity,  therefore  the 
system  is  assumed  to  stay  elastic  throughout  lithiation. 

•  For  galvanostatic  studies,  the  total  current  to  the  Si  NW  was 
maintained  constant;  the  current  density  at  the  surface  of  the  Si 
NW  was  taken  to  be  constant  spatially,  however  it  changes  with 
time  in  accordance  with  the  increase  in  surface  area  related  to 
volume  expansion. 

•  Isothermal  conditions  were  assumed  during  lithiation  of  the  Si 
NW. 


Nus  +  Ns  =  cTv<>  (4) 

The  partial  molar  volume  of  species  LiS,  VLis  is  given  as  a  func¬ 
tion  of  the  molar  volume  of  the  host  material  Vs  and  the  expansion 
factor  f ,  where  ?  can  be  considered  as  an  experimentally  measured 
parameter  defined  as  the  percentage  volume  change  expressed  in 
fraction,  at  a  particular  value  of  A z. 

Vus  =  Vs(l+^)  (5) 


2.2.  Model  equations 
2.2 A.  Material  balance 

The  electrochemical  equilibrium  reaction  between  Li  and  Si  is 
written  in  the  form 

AzLi+  +  A  ze~  +  Si  LiAzSi  i  ( 1 ) 

The  above  reaction  can  be  thought  of  a  single  electron  transfer 
reaction  (Li+  +  e-^Li)  ,  followed  by  lithium  alloying  with  Sii/Ax, 
where  Az  (varies  between  0  and  1 )  is  the  insertion  coefficient  of  Li 
in  LiAzSii/Ax  and  Ax  is  the  maximum  number  of  moles  of  Li  that  can 
reversibly  alloy  with  Si.  The  value  of  Ax  was  measured  to  be  22/5  in 
earlier  work  [32]  at  high  temperatures,  however  at  room  temper¬ 
atures  with  electrochemical  lithiation,  the  value  of  Ax  was 
measured  to  be  15/4  by  different  groups  [33-35].  In  our  work,  the 
value  for  Ax  was  taken  to  be  15/4.  Based  on  the  earlier  work  26], 
the  binary  species  chosen  are  the  empty  (Li  free)  host  lattice  and 
the  lithiated  host  lattice.  Note,  the  host  lattice  in  this  work  is  1  / AxSi 
and  the  lithiated  host  is  LiSii/Ax  and  will  be  denoted  hereafter  as  S 
and  LiS  respectively. 

In  the  Si  NW  electrode,  the  flux  of  LiS  is  obtained  through  the 
generalized  Maxwell-Stefan  equation  [36].  Considering  ideal  so¬ 
lution,  ignoring  thermal  and  external  forced  diffusion  effects  in  the 
generalized  Maxwell-Stefan  relation  and  rearranging  for  the  flux 
of  species  Nus  we  obtain 


Nus  =  *Lis(NuS  +Ns)  -  cTD LiS,S 


aUS^xUS  + 


*LiS 

RT 


(2) 


where  Nus,  Ns  and  Xus,  *s  are  the  molar  fluxes  and  mole  fractions  of 
the  respective  species,  Mus  and  VLis  are  the  molar  mass  and  molar 
volume  of  LiS,  Cj  is  the  total  concentration,  i.e.  Cj  =  cus  +  Cs,  aus  is 
the  thermodynamic  factor  and  p  is  the  density  of  the  material 

p  =  Cj^XiMi  (3) 


where  Az  is  the  insertion  coefficient  of  Li  in  Li^Sii/^.  The  mass 
balance  for  the  species  LiS  is  written  as 


9(ctXus)  +  v.Nljs  =  o  (6) 

The  total  concentration,  cj  and  the  pressure  p  will  be  described 
after  discussion  of  the  diffusion  induced  stresses. 

At  the  outer  boundary  of  the  Si  NW,  a  constant  current  flux 
condition  was  used  as  the  boundary  condition,  while  at  the  center 
an  axial  symmetry  condition  was  used.  In  a  strict  sense,  the  elec¬ 
trode  kinetics  and  the  mass  transfer  of  the  species  in  solution  could 
determine  the  actual  current  distribution  at  the  electrode  surface, 
however  in  this  study,  we  have  restricted  our  simulations  for  the 
case  of  uniform  current  distribution  along  the  nanowire. 


ft'Nuslr=i?NW(t),z 


'ft'N  Lislr,z=HNw(t) 


— SjiappM 
F 

— Sjiapp(t) 
F 


(7) 


~n'Nus\r=0,z  =  0  (8) 

where  iapp(t)  is  the  current  normalized  to  the  external  surface  area 
of  the  Si  NW. 


2.2.2.  Force  balance 

The  strain  for  the  case  of  small  deformation  is  generally 
described  in  the  tensor  notation  as 

€  =  i  (uvx  +  (UVX)T)  (9) 

For  large  deformation  analysis,  the  Lagrangian  strain  (Green 
strain)  is  related  to  the  displacement  gradient  as  follows: 

e  =  i  (uvx  +  Vxu  -  UVX(UVX)T)  (10) 


The  flux  of  S  is  considered  purely  to  be  convective,  and  the 
lattice  velocity  is  defined  through  the  molar  averaged  velocity,  i.e. 
v°  =  xusvus  +  *sVs-  The  total  molar  flux  is  related  to  the  molar 
average  velocity  as 


where  u  is  the  displacement  vector  calculated  from  the  current  and 
the  initial  configuration  of  the  volume  element.  For  finite  defor¬ 
mation,  the  strain  tensor  for  the  2-d  axial  symmetry  case,  in  the 
cylindrical  co-ordinates  is  written  as 


err 
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where  u,  v,  w  are  the  displacements  corresponding  to  elastic  strains, 
in  the  r,  6,  and  z  directions  respectively.  The  displacement  v  in  the  6 
direction  is  zero  based  on  the  axisymmetry  assumption;  subse¬ 
quently  ere  and  eoz  are  also  zero.  The  symmetric  stress  tensor  de¬ 
scribes  the  stress  components  in  the  material  and  contains  three 
normal  stress  Gmaoo,Gzz  and  the  three  symmetric  shear  stresses,  i.e. 
irdJdzJrz  and  the  components  are  given  as 


P  =  ~tr(T)  (19) 

Finally,  the  equilibrium  force  balance  equation  in  the  2- 
d  axisymmetry  geometry  is  expressed  as 


darr  9lrz  (Trr  ~  ^dd 

dr  +  0z  +  r 


(20) 


Grr 

Tr0 

T rz" 

*zr 

Tdr 

?0z 

\  Trd  = 

Tdr 

(12) 

Jzr 

TZ0 

Ozz_ 

{  *rz  = 

'Czr 

The  elastic  stresses  are  correlated  to  the  strains  using  the  elastic 
moduli  matrix,  which  is  a  fourth  order  tensor,  because  of  the 
symmetry  and  isotropic  assumption,  the  number  of  independent 
parameters  in  this  matrix  is  reduced  to  2,  i.e.  the  Lame  parameters, 
A  and  n-  The  stress-strain  relationship  therefore  reduces  to 

x  =  Atr  (e)(1)  +  2/r(e)  (13) 

The  Lame  parameters  could  be  related  to  the  more  commonly 
used  material  properties,  Young’s  modulus  (E)  and  Poisson’s  ratio 
(v)  through  the  relations. 


(3A  +  2 fi)n  (  _  X 

A  +  /i  ’  V  ~  2  (A  +  jii) 


(14) 


In  this  system,  the  insertion  of  Li  into  the  Si  host  introduces  a 
significant  volume  change,  atypical  of  common  insertion  electrodes 
such  as  LiMn204,  LiCoC^,  LiTisO^  etc.,  where  the  total  volume 
change  is  typically  less  than  10%  and  therefore  ignored  in  most 
models.  To  include  the  volume  change  formalism  into  the  modeling 
framework,  the  total  strain  in  the  electrode  is  defined  as  the  sum¬ 
mation  of  the  two  components,  chemical  strain  (stress  free)  and 
elastic  strain,  i.e.,  ej  =  eCh  +  eei;  and  the  chemical  strain  is  expressed 
as  a  function  of  the  partial  molar  volumes  and  mole  fractions  of  the 
species  LiS  and  S,  (note  Xus  +  xs  =  1): 


(Yu s 


(15) 


Consequently,  the  elastic  strain  can  be  written  as  the  difference 
between  the  total  strain  and  the  chemical  strain,  which  can  be 
substituted  back  into  equation  (12)  to  obtain  the  stress— strain 
relationship  for  electrodes  undergoing  elastic  and  chemical  strains. 


e  =  eT  -  ech  (16) 

The  total  concentration  eras  described  in  equation  (6)  is  written 
as  a  function  of  composition  and  pressure,  which  is  related  to  the 
trace  of  the  stress  tensor 


Ct  =  ?(xus,xs)^(tr(x))  (17) 

The  function  f  is  purely  composition  dependent  and  can  be 
defined  as  a  function  of  partial  molar  volumes  of  the  individual 
species.  The  function  \p  is  evaluated  similar  to  the  work  of  Chris¬ 
tensen  et  al.  [27]  using  the  definition  of  a  compressibility  factor  in 
terms  of  differential  volume  element  and  mean  normal  pressure. 
Subsequently,  equation  (17)  is  expressed  as 


Ct  = 


(xUS^US+Xs^s)  lexP 


-3Kfr(t) 


(18) 


where  I<  is  the  bulk  modulus  of  the  material.  The  pressure  defined 
in  equation  (2)  is  the  thermodynamic  pressure,  which  is  assumed 
to  be  equivalent  to  the  mean  normal  pressure,  and  is  expressed  as 


dTrz  dOzz  Trz  =  Q 

dr  +  9z  +  r 


(21) 


The  individual  velocity  components  of  v°  in  equation  (4)  is 
calculated  from  the  time  derivatives  of  corresponding  displace¬ 
ment  fields. 


r,z 


(22) 


aw 
Vz  =  W 


(23) 


Equations  (6),  (18),  (20)  and  (21)  were  used  to  solve  for  the 
variables  u,  w,  xus,  and  ct  respectively  in  the  Si  NW  domain.  The 
equilibrium  force  balance  equations  (20)  and  (21)  are  also  valid  in 
the  Cu  CC  domain,  with  the  exception  that  the  total  strain,  ej  is 
purely  elastic  and  the  chemical  strain  component  is  absent  based 
on  the  assumption  that  lithium  does  not  transport  into  the  Cu 
substrate.  Therefore  the  displacement  components  u,  w  are  the 
only  variables  to  be  solved  for  in  the  Cu  CC  domain. 

The  force  balance  equation  in  the  Si  NW  was  constrained  to  the 
following  boundary  conditions 


I  azz\ r,z=HNW(t)  =  0 


(24) 


The  base  of  the  Cu  CC  substrate  is  subjected  to  a  fixed  constraint 
boundary  conditions 


f  ulr,z=0  =  0 
\w|r,z=o  =  0 


(25) 


The  outer  surface  of  the  Si  NW,  the  top  unanchored  portion  of 
the  Cu  CC  and  the  outer  surface  of  the  Cu  CC  are  assumed  to  be  free 
surfaces  i.e. 


n-T  =  0 


(26) 


3.  Solving  methodology  &  parameters 

The  equations  are  solved  using  COMSOL  Multiphysics  with  the 
structural  mechanics  module  to  solve  for  the  displacement  com¬ 
ponents  u,  w  equations  (20)  and  (21 )  and  a  general  PDE  interface  to 
solve  for  xys  equation  (6).  The  mass  balance  equations  were  re¬ 
written  in  terms  of  material  derivatives  for  the  ease  of  imple¬ 
mentation  in  the  material  framework  in  COMSOL. 

Since  the  dimensions  of  the  Si  NW  change  significantly  upon 
lithiation,  the  initial  mesh  configuration  has  to  be  updated  at  each 
time  step  to  accommodate  for  the  updated  geometry.  In  this  model, 
the  deformation  of  the  mesh  is  determined  by  the  displacement 
components  (u,  w)  calculated  from  structural  mechanics  module. 
The  technique  for  mesh  movement  is  called  an  Arbitrary 
Lagrangian— Eulerian  (ALE)  method,  which  is  an  intermediate  be¬ 
tween  the  Lagrangian  and  Eulerian  methods,  and  it  allows  moving 
boundaries  without  the  need  for  the  mesh  movement  to  follow  the 
material.  For  the  2-d  model,  the  Si  NW  part  of  the  geometry  was 
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Table  1 

List  of  parameter  values  used  in  the  simulation. 


Parameter 

Value 

Units 

Partial  molar  volume  of  LiS,  VLis 

13.11 

ml  mol”1 

Molar  volume  of  S,  Vs 

3.214 

ml  mol”1 

Maximum  number  of  moles  of  Li  that  can 

15/4 

No  units 

reversibly  alloy  per  mole  of  Si,  Ax 

Maximum  insertion  coefficient  of  Li  in  LiS,  Azmax 

1 

No  units 

Diffusion  coefficient  of  Li  in  Si,  DLis,s 

2  x  10-12 

cm2  s”1 

Young’s  modulus  of  LiS,  E 

92.16 

GPa 

Poisson’s  ratio  of  LiS,  v 

0.27 

No  units 

Young’s  modulus  of  Cu  substrate,  ECu 

110 

GPa 

Poisson’s  ratio  of  Cu  substrate,  vCu 

0.35 

No  units 

Thermodynamic  factor,  a 

1 

No  units 

Expansion  faction,  f  (measured  at  Azmax) 

3.079 

No  units 

mapped  with  300  mode  points  along  the  axial  direction,  and  100 
points  along  the  radius,  while  the  base  Cu  CC  structure  was  mapped 
with  150  node  points  along  the  axial  direction  and  150  node  points 
along  the  radius.  In  all,  the  geometry  consisted  of  52,500  quadri¬ 
lateral,  1300  line  and  7  vertex  elements.  An  absolute  tolerance  of 
1CT15  and  1CT6  was  used  for  the  displacement  variables  (u,  w)  and 
Xus  respectively,  and  a  relative  tolerance  of  1CT5  was  used  to 
establish  convergence.  Automatic  time  stepping  (based  on  the 
solver)  was  used,  and  the  computational  run  time  taken  for  a 
complete  charge  (167  time  steps)  was  17281  s,  using  a  16  core  Intel 
Xeon  2.27  GHz  processor.  The  parameters  used  in  the  model  are 
given  in  Table  1. 

4.  Results  and  discussion 

4.1  Diffusion  induced  stresses 

Fig.  2a  shows  the  mole  fraction  distribution  of  LiS  in  the  Si  NW  at 
the  end  of  lithiation.  The  lithiation  rate  in  this  simulation 


corresponds  to  a  surface  current  density  of  0.02  mA  cm~2( initial) 
equivalent  to  a  1  -h  rate.  The  solid  line  in  the  plot  marks  the  initial 
undeformed  configuration  of  the  Si  NW  anchored  to  the  Cu  CC 
substrate.  In  this  study,  the  simulation  was  terminated  when  the 
local  mole  fraction  reached,  Xus  =  1  anywhere  within  the  electrode. 
As  observed  from  the  plot,  the  top  surface  of  the  Si  NW  is  mass 
transfer  limited  and  gets  completely  lithiated  (xus  =  1 )  while  the 
bulk  of  the  Si  NW  is  still  partially  lithiated  (xus  =  0.86),  which  limits 
complete  electrode  utilization.  The  final  deformed  configuration  of 
the  Si  NW  shows  the  increase  in  the  radial  and  axial  dimensions  of 
the  Si  NW  due  to  the  combination  of  chemical  and  elastic  strains 
during  lithiation.  The  top  of  the  Si  NW  is  expanded  more,  due  to 
maximum  lithiation  in  that  region  resulting  in  increased  chemical 
strain,  and  regions  very  close  to  the  Si  NW/Cu  CC  interface  (as 
shown  in  Fig.  2b)  are  minimally  expanded  as  the  lithiation  is 
limited  due  to  the  high  stresses  developed  at  the  lithium  blocking 
interface.  Fig.  2b  also  shows  the  displacement  of  the  Si  NW/Cu  CC 
interface  from  the  initial  configuration  due  to  the  interfacial 
stresses.  Note,  the  Si  NW  region  is  pushed  into  the  Cu  CC  region  (z 
axis)  by  - 1  nm. 

Fig.  3  shows  a  snap  shot  of  the  local  volumetric  strain  distri¬ 
bution  in  the  structure  at  the  end  of  lithiation.  The  volumetric  strain 
in  the  Si  NW  is  non-uniform  in  the  axial  direction,  especially  at  the 
top  and  the  Si  NW/Cu  CC  interface  regions.  In  general,  the  local 
volumetric  strain  distribution  in  the  Si  NW  domain  correlates  to  the 
concentration  distribution  in  Fig.  2a,  because  the  total  volumetric 
strain,  defined  in  equation  (16)  is  predominantly  determined  by  the 
chemical  strain.  In  the  Cu  domain,  the  volumetric  strains  are  purely 
elastic  and  mostly  tensile,  except  at  regions  close  to  the  Si  NW/Cu 
CC  and  away  from  the  center  where  some  regions  are  compressive. 

Fig.  4  compares  the  radial,  tangential  and  axial  stress  compo¬ 
nents  across  the  radius  of  the  Si  NW  at  different  times  during 
lithiation.  The  radial  cut  section  in  this  plot  is  taken  at  half  the 
initial  height  of  the  Si  NW  (z  =  Hqu  +  Hnw/2).  Several  features  are 
observed  in  this  plot;  firstly  the  radial  stress  across  the  radius  of  the 
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(a)  Mole  fraction  profile  of  LiS,  i .e.XLis  in  the  entire  Si  NW  (b)  Mole  fraction  profile  of  LiS,  l.e.XLis  close  to 
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Fig.  2.  Mole  fraction  distribution  in  the  Si  NW  at  end  of  lithiation.  A  constant  1-hr  rate  current  equivalent  to  a  current  density  (initial)  of  0.0208  mA  cm  2  was  used  in  this 
simulation. 
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Fig.  3.  Volumetric  strain  distribution  in  the  Si  NW/Cu  CC  at  the  end  of  lithiation.  A 
constant  1-hr  rate  current  equivalent  to  a  current  density  (initial)  of  0.0208  mA  cm-2 
was  used  in  this  simulation. 


Si  NW  is  always  tensile,  and  is  maximum  at  the  center.  This  is 
because  of  the  concentration  build  up  at  the  surface  which  causes 
volumetric  strain  in  the  outer  layers,  which  in  turn  radially  pulls  the 
inner  layers  to  create  the  tensile  stresses  throughout  the  radius  of 
the  Si  NW  (as  plotted  in  Fig.  4a,  b,  c,  &  d).  As  a  function  of  time,  the 
maximum  radial  stresses  (at  the  center  of  the  Si  NW)  increase  up  to 
the  first  10  s  ( ~43  MPa  observed  at  10  s)  and  continuously  decrease 
at  longer  times  (~2  MPa  observed  at  1000  s).  Also,  the  maximum 
tangential  and  the  axial  stresses  follow  a  similar  trend.  This 
behavior  is  due  to  the  competing  effects  of  the  chemical  diffusion 
term  and  the  pressure  induced  term  in  equation  (2)  towards  the 
overall  flux  of  the  species.  At  short  times,  the  species  flux  is 
dominated  by  the  chemical  diffusion  term,  while  at  longer  times,  is 
dominated  by  pressure  gradient  term,  resulting  in  reduced  con¬ 
centration  gradients.  Since  the  build  up  of  stresses  is  proportional 
to  concentration  gradients  based  on  equation  (16)  and  the  force 
balance  relations  (equations  (20)  and  (21)),  the  stresses  decrease 
when  the  flux  is  dominated  by  the  pressure  gradients.  Secondly,  the 
tangential  and  the  axial  stresses  are  always  compressive  towards 
the  outer  surface  and  tensile  towards  the  inner  core.  This  behavior 
is  due  to  the  radial  expansion  of  the  outer  surface  which  creates 
compressive  strains  in  the  tangential  and  the  axial  direction  to¬ 
wards  the  outer  surface,  while  the  inner  core  is  pulled  outwards  in 
all  coordinates  creating  tensile  stresses  in  the  tangential  and  axial 
directions.  Thirdly,  at  the  center  of  the  Si  NW,  the  radial  and  the 
tangential  stresses  are  equal,  and  the  magnitude  of  the  axial  stress 
is  twice  the  radial  or  tangential  stress.  This  scenario  is  represen¬ 
tative  of  a  1  -d  plane  strain  condition  with  infinitesimal  deforma¬ 
tion,  where  the  condition  an  =  aod  =  &zzl 2  is  satisfied  at  the  center 
(r  =  0)  in  cylindrical  coordinates.  This  behavior  suggests  that  far 


Fig.  4.  Profiles  of  radial,  tangential,  and  the  axial  stresses  across  the  radius  of  the  Si  NW  at  half  the  height  of  the  nanowire  (z  =  HCu  +  HNW/ 2)  at  various  times:  a)  1  s,  b)  10  s,  c)  100  s 
&  d)  1000  s.  A  constant  1-hr  rate  current  equivalent  to  a  current  density  (initial)  of  0.0208  mA  cm-2  was  used  in  these  simulations. 
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Fig.  5.  Profiles  of  radial,  tangential,  and  the  axial  stresses  across  the  radius  of  the  Si  NW,  close  to  the  Si  NW/Cu  CC  interface  (z  =  HCu)  at  various  times:  a)  1  s,  b)  10  s,  c)  100  s  &  d) 
1000  s.  A  constant  1-hr  rate  current  equivalent  to  a  current  density  (initial)  of  0.0208  mA  cm-2  was  used  in  these  simulations. 


away  from  the  Si  NW/Cu  CC  interface,  the  stress  behavior  is  similar 
to  a  1-d  plane  strain  condition. 

Fig.  5  compares  the  radial,  tangential  and  axial  stress  compo¬ 
nents  close  to  the  Si  NW/Cu  CC  interface  on  the  Si  NW  side.  At 


Fig.  6.  Evolution  of  maxim  radial  stress  (r  =  0),  maximum  tangential  stress  (r  =  i?Nw(t)) 
and  maximum  axial  stress  (r  =  0)  with  time  at  half  the  height  of  the  nanowire 
(z  =  Hc u  +  Hnw/2).  A  constant  1-hr  rate  current  equivalent  to  a  current  density  (initial) 
of  0.0208  mA  cm-2  was  used  in  this  simulation.  Plot  in  the  subset  shows  the  evolution 
of  the  ratio  of  diffusive  flux  and  the  pressure  induced  flux  to  the  overall  surface  flux 
with  time.  A  constant  1-hr  rate  current  equivalent  to  a  current  density  (initial)  of 
0.0208  mA  cm-2  was  used  in  these  simulations. 


short  times,  t  =  1  and  t  =  10  s,  the  stress  profiles  across  the  radius 
matched  quantitatively  with  the  stresses  profiles  across  the  radius 
in  the  center  region  of  the  Si  NW(as  discussed  in  the  earlier  sec¬ 
tion),  while  at  longer  times  the  presence  of  the  constraint  (sub¬ 
strate)  significantly  alters  the  stress  profiles.  At  100  s,  the  radial 
and  tangential  stress  profiles  are  similar  to  short  time  behavior, 
however  the  axial  stress  becomes  less  tensile  at  the  center  and  at 
longer  times  (t  =  1000  s)  the  axial  stress  reverses  its  general  trend 
and  becomes  compressive  in  the  inner  part  and  tensile  at  the 
outer  part.  Furthermore,  the  magnitude  of  all  the  stress  compo¬ 
nents  increases  considerably  and  are  in  the  range  of  250- 
500  MPa.  This  is  possibly  because  of  the  constant  lithiation  flux 
(imposed  by  the  boundary  conditions)  in  the  regions  close  to  the 
interface,  while  simultaneously  the  interface  is  also  being  con¬ 
strained  by  the  Cu  substrate  resulting  in  the  large  tensile  and 
compressive  stress  regions.  These  enormous  stress  components 
could  potentially  cause  the  Si  NW  structure  to  yield  or  fracture  in 
these  regions  close  to  the  interface. 

Fig.  6  shows  the  evolution  of  maximum  stress  for  each 
component  with  time  during  lithiation  of  the  Si  NW.  In  this  plot,  the 
z  co-ordinate  is  halfway  through  the  initial  height  of  Si  NW,  i.e. 
z  =  He u  +  Hnw/2  and  the  r  co-ordinate  is  chosen  corresponding  to 
where  the  maximum  value  of  stress  in  each  component  occurs.  The 
maximum  radial  and  axial  stresses  are  tensile  and  always  occur  at 
the  center  (r  =  0)  of  the  Si  NW,  while  the  maximum  tangential 
stresses  are  compressive  and  occur  at  the  outer  surface  of  the  Si 
NW,  i.e.  r  =  f?Nw(t)-  For  the  parameters  used  in  the  simulation,  the 
values  for  the  stress  components  peak  at  -4  s  and  decrease  sub¬ 
sequently.  As  explained  in  the  earlier  section,  this  behavior  is  due  to 
the  competing  effects  of  the  chemical  diffusion  term  and  the 
pressure  induced  term  in  equation  (2)  towards  the  overall  flux  of 
the  species.  The  plot  in  the  subset  of  Fig.  6  clearly  shows  the  shift 
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(a)  Radial  stresses  at  the  center  of  the  Si  NW,  i.e  at  (z,r)  = 

(HCu  +  UfSL,0) 


Time  (s) 


(b)  Tangential  stresses  at  the  outer  radius  of  the  Si  NW,  i.e.  at  (z,r)  = 

(■ HCu  + 

Fig.  7.  Evolution  of  stresses  for  different  lithiation  rates.  1C  rate  corresponds  to  a 
constant  1-hr  rate  current  equivalent  to  a  current  density  (initial)  of  0.0208  mA  cm-2. 

from  the  diffusion  dominated  transport  at  short  times,  to  pressure 
driven  transport  at  longer  times. 

4.2.  Effect  of  lithiation  rate 

Figs.  7a  and  b  shows  the  effect  of  lithiation  rate  on  the  evolution 
of  the  maximum  radial  and  tangential  stresses  with  time.  Here  C 
rate  corresponds  to  a  1-h  rate  equivalent  to  an  initial  current 
density  of  0.02  mA  cm-2.  The  inset  plots  in  Figs.  7a  and  b  show  a 
linear  increase  in  maximum  radial  and  tangential  stresses  with 
lithiation  rate.  Also  at  higher  lithiation  rates,  the  peak  maximum 
stresses  occur  at  shorter  times  as  seen  from  the  shift  in  the  peak 
towards  the  left.  This  behavior  suggests  possibility  for  mechanical 
fracture  at  very  short  times  under  high  current  conditions,  typically 
seen  in  hybrid  electric  vehicle,  fast  charge  or  regenerative  braking 
applications,  despite  the  nanoscale  dimensions  of  the  electrode. 
The  occurrence  of  peak  maximum  stresses  at  shorter  times  at 


higher  lithiation  rates,  could  also  be  explained  through  the  inter¬ 
play  between  the  diffusive  and  the  pressure  induced  flux.  At  higher 
rates,  large  concentration  gradients  (due  to  chemical  diffusion)  are 
established  at  shorter  times,  creating  a  large  stress  field  at  the 
surface.  Consequently,  a  larger  pressure  gradient  is  built  up,  which 
then  dominates  the  species  flux,  compared  to  the  chemical  diffu¬ 
sion  mode  in  the  bulk  of  the  Si  NW. 

4.3.  Effect  of  Si  NW  radius 

Figs.  8a  and  b  show  the  effect  of  Si  NW  radius  on  the  evolution  of 
maximum  radial  and  tangential  stresses  with  time  for  1-h  lithiation 
rate.  The  1-h  lithiation  rate  (C  rate)  for  the  Si  NW  radii  of  50, 100, 
200,  300  and  500  nm  corresponds  to  a  current  density  of  0.0208, 
0.0415  and  0.0826,  0.0826,  and  0.1251  mA  cm-2  respectively.  The 


(a)  Radial  stresses  at  the  center  of  the  Si  NW  i.e.  at  (z,r)  = 

{Hou  +  k f^,0) 


(b)  Tangential  stresses  at  the  outer  radius  of  the  Si  NW,  i.e.  at  (z,  r)  = 
(■ HCu  +  ,  Rnw  (t)) 


Fig.  8.  Evolution  of  stresses  for  different  radii  of  Si  NW  at  1-hr  lithiation  rate.  The 
current  used  in  each  case  (C  rate),  corresponds  to  a  current  density  (initial)  equivalent 
to  0.0208  mA  cm-2,  0.0415  mA  cm-2,  0.0826  mA  cm"2,  0.1251  mA  cm-2,  and 
0.2085  mA  cm-2,  for  Si  NW  of  radius  50, 100,  200,  300  and  500  nm  respectively. 
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increased  surface  current  density  explains  the  higher  radial  and 
tangential  stresses  observed  in  structures  with  larger  radii.  Also, 
the  increased  current  densities  (for  larger  radius),  shift  the  peak 
maximum  stresses  to  longer  times,  which  is  contrary  to  the  effect 
observed  at  higher  current  densities  for  constant  radius  (Fig.  7a  and 
b).  This  behavior  suggests  that  the  increase  in  the  current  density 
(for  larger  radius  structures)  is  not  large  enough  to  counter  the 
longer  diffusion  length,  which  in  turn  delays  the  time  for  maximum 
stresses  to  develop.  Consequently,  the  contribution  from  pressure 
induced  flux  takes  a  longer  time  to  offset  the  diffusion  dominated 
flux  for  nanowires  with  larger  radii.  Further,  since  the  stress  values 
for  particles  of  larger  radius  are  at  any  time  higher  than  that  for  the 
particles  with  smaller  radius,  the  latter  is  preferred  especially  for 
high  rate  applications.  While  smaller  particles  are  not  preferred  due 
to  lower  compressed  density  and  higher  exposed  surface  area  to  the 
electrolyte,  they  clearly  offer  an  advantage  from  a  mechanical  stand 
point.  Design  of  optimal  particle  size  should  however  be  considered 
based  on  the  energy  and  power  requirements  for  specific 
applications. 

5.  Conclusions 

A  2-d  transient  numerical  model  to  simulate  the  electro¬ 
chemical  lithiation  of  Si  NW  is  presented.  The  model  predicts  non- 
uniform  volumetric  strain  along  the  length  of  the  Si  NW,  with  re¬ 
gions  of  maximum  expansion  at  the  top  of  the  Si  NW  and  almost  no 
expansion  close  to  the  Cu  CC  interface.  The  magnitudes  of  the  stress 
components  are  very  high  at  the  Si  NW/Cu  CC  interface,  compared 
to  the  stresses  developed  far  away  from  the  interface.  The  stress 
evolution  with  time  is  strongly  dependent  on  the  relative  magni¬ 
tudes  of  chemical  and  the  pressure  diffusion  fluxes.  The  maximum 
stresses  occur  during  the  time  when  the  flux  is  dominated  by  the 
chemical  diffusion  term,  i.e.  ~1— 10  s  for  the  rates  and  radius 
chosen  for  the  simulations.  Increase  in  radius  of  the  nanowire  and 
increase  in  lithiation  rates  develop  larger  radial  and  tangential 
stresses.  Further,  the  peak  maximum  stresses  occur  at  shorter  times 
with  increase  in  lithiation  rates,  while  it  occur  at  longer  times  with 
increase  in  the  radius  of  the  nanowire.  The  comparison  of  the  re¬ 
sults  from  the  2-d  numerical  model  to  the  1-d  (radial)  model  and 
the  evaluation  of  stresses  for  different  geometries  are  presented  in 
Part  II  of  this  work. 
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Nomenclature 

c  concentration  of  species,  Cs,  cus  or  ct,  mol  m-3 

Nl is  flux  of  species  Cs,  or  cy s  mol  m”2  s-1 

x  mole  fraction  of  species,  Xs  or  xys 

Az  insertion  coefficient  of  Li  in  LiAzSii/A* 

Ax  maximum  number  of  moles  of  Li  that  can  reversibly  alloy 
per  mole  of  Si 
t  time,  s 

r  radial  co-ordinate,  m 

6  tangential  co-ordinate,  deg 

z  axial  co-ordinates,  m 

v  velocity  of  species,  vs  or  vys,  m  s-1 

v°  molar  average  velocity,  m  s-1 

Dys,s  binary  diffusion  coefficient  of  Li  in  LiS,  m2  s-1 
aus  thermodynamic  factor 

M  molar  mass  of  species,  Ms  or  Mys,  g  mol-1 

F  Faraday’s  constant,  C  g-1  equiv 

R  universal  gas  constant,  J  mor1  K1 


T  temperature,  K 

Si  stoichiometric  coefficient 

V  partial  molar  volume  of  species  LiS,\/Lis  or  molar  volume 

of  host  material  S,  Vs ,  m3  mol-1 
?  expansion  factor  as  described  in  equation  (5) 

p  density  of  LizSii/A*,  g  m~3 

p  thermodynamic  pressure,  N  m-2 

u  displacement  vector  field 

u  displacement  in  r  co-ordinate,  m 

v  displacement  in  6  co-ordinate,  m 

w  displacement  in  z  co-ordinate,  m 

n  unit  outward  normal  vector,  m 

z'app  current  density  related  to  the  outer  surface  area  of  the  Si 
NW,  A  m”2 

Fnw  radius  of  the  nanowire,  also  a  function  axial  co-ordinate 
and  time,  m 

H  height  (z  axis)  of  the  structure,  i.e.  nanowire  or  Cu  current 

collector 

t  stress  tensor  matrix 

a  stress  components  of  the  stress  tensor  matrix,  N  m-2 

e  strain  tensor  matrix 

6  strain  components  of  the  strain  tensor  matrix,  N  m-2 

X  Lame  constant  given  by  Ev\\  +  v,  N  m-2 

p  Lame  constant  given  by  E/2(l  +  r)(l  -  2v ),  N  m-2 

E  Young’s  modulus,  N  m-2 

I<  isothermal  bulk  modulus  given  by  F/ 3(1  -  2v ),  N  m-2 

v  Poisson’s  ratio 

Subscripts 

LiS  lithiated  host 

S  unoccupied  host 

ch  chemical 

T  total 

max  maximum 

NW  nanowire 

Cu  Cu  current  collector 
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